#ifndef PBVP_DUAL_NODE_H
#define PBVP_DUAL_NODE_H

#include <vector>
#include <iostream>
#include <algorithm>
#include <QString>

#include "moment.h"
#include "attribute.h"

#include "../contree/augmentedctree.h"

namespace PBVP {
	class GeometryData;

	/// Similarity score data structure
	/**
	* It has a separate field for each set of attributes
	*/
	class SimScore {
	public:
		/// global weights for similarity score computation
		static float wGeo, wShape, wPot, wTopo;

		/// four channel of scores
		float geoScore, shapeScore, potScore, topoScore;

		/// Default constructor
		SimScore () {
			geoScore = shapeScore = potScore = topoScore = 0;
		}

		/// multiply a constant
		SimScore operator * (float f) const {
			SimScore sim;
			sim.geoScore = geoScore * f;
			sim.shapeScore = shapeScore * f;
			sim.potScore = potScore * f;
			sim.topoScore = topoScore * f;
			return sim;
		}   

		/// divide a constant
		SimScore operator / (float f) const {
			assert(f > 0);
			SimScore sim;
			sim.geoScore = geoScore / f;
			sim.shapeScore = shapeScore / f;
			sim.potScore = potScore / f;
			sim.topoScore = topoScore / f;
			return sim;
		}

		/// += operator
		const SimScore& operator += (const SimScore& sim) 
		{
			geoScore += sim.geoScore;
			shapeScore += sim.shapeScore;
			potScore += sim.potScore;
			topoScore += sim.topoScore;
			return *this;
		}

		/// weighted average of the individual scores
		float average() const
		{
			return(wGeo*geoScore + wShape*shapeScore + wTopo*topoScore + wPot*potScore);
		}

		/// Set the weights of the similarities of individual attributes
		/**
		* @note static function
		*/
		static void setWeights(float _wGeo, float _wShape, float _wTopo, float _wPot)
		{
			wGeo = _wGeo;
			wShape = _wShape;
			wTopo = _wTopo;
			wPot = _wPot;
		}

		/// Print the individual simscores
		void print() 
		{
			printf("geoScore = %f, shapeScore = %f, topoScore = %f, potScore = %f, average score = %f\n",
				   geoScore, shapeScore, topoScore, potScore, average());
		}

		/// Convert to QString
		QString toQString() const;

	};

	/// Dual Contour Tree (DCT) Node
	/**
	*	It has a vector of attributes that can be used for comparison
	*/
	class DualNode {
	public:
		/// unique ID of the node
		int id;
		/// functional range 
		int range;
		/// number of critical points in the interval volume
		int ncrit;  
		// max and min func values of this node
		float fmax, fmin;
		/// normalized volume corresponding to this node;
		float normVolume;   
		/// unnormalized (raw) volume
		float rawVolume;            

		///
		float fint;							// Integral of a second function 
		BettiNumber bl, bu;					// lower and upper boundary Betti numbers

		/// lower neighbors of this node
		std::vector<DualNode *> *p_below;
		/// upper neighbors of this node
		std::vector<DualNode *> *p_above;
		DualNode* parent;					// points to the dual node in coarser level
		DualNode* match;					// A dual node in another graph that matches this node.		
		/// vector of attributes									 
		std::vector<Attribute *> attributes;

		/// Id of the tetra in the DCT node as start for marching
		int startTetra;

		///Default Constructor
		/**
		* @param _range: interval range of the node
		*/
		DualNode(int _range = 0); 

		/// Destructor
		~DualNode(); 

		/// Compute the attributes of the DCT node
		/**
		* 
		* @param data: 3D regular volumetric data
		* @param start_tid: the id of the first tetrahedron from which
		*           marching starts.
		*/
		// void computeAttributes(VolumeReg3Critical* data, int start_tid);

		/// Extract bounding surfaces
		// Surface3D* extractSurface(VolumeReg3Critical* data, bool upperSurface = true);

		/// add an attribute to the attribute vector
		void addAttribute(Attribute* attr)
		{
			attributes.push_back(attr);
		}

		/// Get the volume of this node that is normalized by the total volume of the DCT
		float getNormalizedVolume() const {
			return normVolume;
		}

		/// Get the volume of this node
		float getVolume() {
			GeoAttribute* pGeo = (GeoAttribute *)attributes[0];
			return pGeo->volume;
			//return rawVolume;
		}

		/// Get the Integral of a property function
		float getFuncIntegral(int ip = 0);

		/**
		* Merge a lower level dual node to a higher level node
		* @param direct: 0 -- merge from under, 1 -- merge from above
		*/
		void merge(DualNode* p_node, int direct,  bool firstChild = false);

		///	Add a neighbor of the lower range
		void addLowerNeighbor(DualNode* node) {
			p_below->push_back(node);
		}

		/// Add a neighbor of the upper range.
		void addUpperNeighbor(DualNode* node) {
			p_above->push_back(node);
		}

		/// The number of neighbors in the lower range.
		int lowerDegree() const {
			return(int) (p_below->size());
		}

		///	The number of neighbors in the upper range.
		int upperDegree() const {
			return int (p_above->size());
		}

		/// Get the ith neighbor in the lower range.
		DualNode * lowerNeighbor(int i) const {
			return(*p_below)[i];
		}

		/// Get the ith neighbor in the lower range.
		DualNode * lowerNeighbor(int i) {
			return(*p_below)[i];
		}

		/// Get the ith neighbor in the upper range.
		DualNode * upperNeighbor(int i) const {
			return(*p_above)[i];
		}

		/// Get the ith neighbor in the upper range.
		DualNode * upperNeighbor(int i) {
			return(*p_above)[i];
		}

		/// Check if a node has additional property attributes
		bool hasPropAttributes() const {
			return attributes.size() > 3;
		}

		/// Get the number of additional property attributes
		int numOfPropAttributes() const {
			return attributes.size() - 3;
		}


		///	Get the importance of this node 
		float importance() {
			//return 0.5*(fmax-fmin) + 0.5*volume;
			return 0.1f*(fmax-fmin) + 0.9f*normVolume;
		}

		/// Set Max functional value of the DCT node
		void setFuncMax(float x) {
			if (x > fmax) fmax = x;
		}

		/// Set Min functional value of the DCT node
		void setFuncMin(float x) {
			if (x < fmin) fmin = x;
		}

		///	Matchup this node with another node
		void setMatchNode(DualNode* pNode) {
			match = pNode;
			pNode->match = this;
		}

		///  Check if this node has a match node
		bool isMatched() {
			return(match != NULL);
		}

		/// clear the matched node
		void clearMatch() {
			match = NULL;
		}

		/**
		* Check if the parent of this node has a match node
		*/
		bool isParentMatched(DualNode* pNode) {
			if (parent == NULL)	return true; // top level node
			//return (parent == NULL || parent->match != NULL);
			return(parent->match != NULL && parent->match == pNode->parent);
		}

		/// Check if the node has it attributes computed
		bool isVolumeComputed() const {
			return volumeComputed;
			// return(normVolume >= 0);
		}

		/// Set the attribute computed flag to true
		void setVolumeComputed() {
			volumeComputed = true;
		}

		/// Connect a lower noder with an upper node
		static void connect(DualNode* node1, DualNode *node2) {
			node1->addUpperNeighbor(node2);
			node2->addLowerNeighbor(node1);
		}

		/// Compute the individual similarity scores with the attribute vector
		SimScore similarity(DualNode* pNode);

		/// ostream output operator for DualNode.
		friend std::ostream& operator << (std::ostream& out, const DualNode& node) {
			int k;
			out << node.id << " " << node.range << std::endl;
			out << node.fmin << " " << node.fmax << " " << node.rawVolume << " " << node.ncrit << " " << node.fint << std::endl;
			out << node.bl << node.bu << std::endl;
			// out << node.moments;
			// export the attribute vector
			// Four attributes: geometry, shape, property, and topology
			std::vector<Attribute*>::const_iterator it = node.attributes.begin();
			// Attribute 1: geometry
			GeoAttribute* pGeo = (GeoAttribute*)(*it);
			out << *pGeo;
			++it;
			// Attribute 2: shape
			ShapeAttribute* pShape = (ShapeAttribute*)(*it);
			out << *pShape;
			++it;
			// Attribute 3: topology
			TopoAttribute* pTopo = (TopoAttribute*)(*it);
			out << *pTopo;
			++it;

			// Attribute 4: potentials (there may be multiple potentials)
			int nPot = node.attributes.size() - 3;
			out << nPot << std::endl;
			for (int i = 0; i < nPot; i++) {
				PotentialAttribute* pPot = (PotentialAttribute*)(*it);
				out << *pPot;
				++it;
			}

			out << node.startTetra << std::endl;
			out << node.lowerDegree() << " " << node.upperDegree() << " ";
			for (k = 0; k < node.lowerDegree(); k++) {
				out << (node.lowerNeighbor(k))->id << " ";
			}
			for (k = 0; k < node.upperDegree(); k++) {
				out << node.upperNeighbor(k)->id << " ";
			}
			out << std::endl << std::endl;
			return out;
		}


		/// istream input opeartor for DualNode
		friend std::istream& operator >> (std::istream& in, DualNode& node) {
			in >> node.id >> node.range;
			in >> node.fmin >> node.fmax >> node.rawVolume >> node.ncrit >> node.fint;
			in >> node.bl >> node.bu;
			// in >> node.moments;

			// import the attribute vector 
			// Attributes: geometry, shape, topology, and properties 
			// Attribute 1: geometry
			GeoAttribute* pGeo = new GeoAttribute;
			in >> *pGeo;
			node.addAttribute(pGeo);
			// Attribute 2: shape
			ShapeAttribute* pShape = new ShapeAttribute;
			in >> *pShape;
			node.addAttribute(pShape);
			// Attribute 3: topology
			TopoAttribute* pTopo = new TopoAttribute;
			in >> *pTopo;
			node.addAttribute(pTopo);
			// Attribute 4: Potentials
			int nPot;
			in >> nPot;
			for (int i = 0; i < nPot; i++) {
				PotentialAttribute* pPot = new PotentialAttribute;
				in >> *pPot;
				node.addAttribute(pPot);
			}

			in >> node.startTetra;
			return in;
		} 


		/// Print the info about DualNode
		void print() const; 

		/// Convert node info to QString
		QString toQString() const;

		/// Simple score used to pick matched nodes
		float simpleScore(DualNode* pNode) {
			float vsim = 0, fisim = 0;
			if (std::max(normVolume, pNode->normVolume) == 0) {
				vsim = 0;
			} else {
				vsim = std::min(normVolume, pNode->normVolume) / std::max(normVolume, pNode->normVolume);
			}
			if (hasPropAttributes() && pNode->hasPropAttributes()) {
				float f1 = getFuncIntegral() / getVolume(); 
				float f2 = pNode->getFuncIntegral() / pNode->getVolume(); 
				float len = std::max(1.0f, (float)std::max(fabs(f1), fabs(f2)));
				fisim = 1 - fabs(f2-f1)/len;
			}

			return 0.5f*vsim + 0.5f*fisim;
		}

		/// Used when DCT prune	
		bool pruned;

	private:
		/// global min and max
		static float all_min, all_max;      
		/// 
		bool volumeComputed;    
	};

	class Less_DualNodePtr {
	public:
		bool operator() (DualNode *pNode1, DualNode *pNode2)
		{
			return(pNode1->getVolume() < pNode2->getVolume());
		}
	};
}

#endif
